% Marshallian and Hicksian
 axoptions={'scaled ticks = false',...
           'y tick label style={/pgf/number format/.cd, fixed, fixed zerofill,precision=3, set thousands separator={}}',...
           'x tick label style={/pgf/number format/.cd, set thousands separator={}}',...
           'legend style={font=\normalsize}'}; %\scriptsize
 
U_0 = U_vec(:,:,1);
U_T = U_vec(:,:,T);
I_T = I_vec(:,:,T);
I_T(isnan(U_T)) = [];
U_T(isnan(U_T)) = [];
U_0(isnan(U_0)) = [];
UT90 =quantile(U_T,0.90);
UTMed =quantile(U_T,0.5);
UT10 =quantile(U_T,0.1);
U_0mean = mean(U_0);

F_UI_T = interpNaN(log(U_vec(:,:,T)),log(I_vec(:,:,T)));

IT90 = exp(F_UI_T(log(UT90))); 
ITMed = exp(F_UI_T(log(UTMed))); 
IT10 = exp(F_UI_T(log(UT10))); 

%%
bt_U_med = zeros(1,T);
bt_U_10 = zeros(1,T);
bt_U_90 = zeros(1,T);
bt_U_mean0  = zeros(1,T);

bt_I_med = zeros(1,T);
bt_I_10 = zeros(1,T);
bt_I_90 = zeros(1,T);
bt_I_mean0  = zeros(1,T);

for t = 1:T
Fbt_U = interpNaN(log(U_vec(:,:,t)),Bx_vec(:,:,t));
bt_U_med(:,t) = Fbt_U(log(UTMed));
bt_U_10(:,t) = Fbt_U(log(UT10));
bt_U_90(:,t) = Fbt_U(log(UT90));
bt_U_mean0(:,t) = Fbt_U(log(U_0mean));
end

for t = 1:T
Fbt_I = interpNaN(log(I_vec(:,:,t)),Bx_vec(:,:,t));
bt_I_med(:,t) = Fbt_I(log(ITMed));
bt_I_10(:,t) = Fbt_I(log(IT10));
bt_I_90(:,t) = Fbt_I(log(IT90));
bt_I_mean0(:,t) = Fbt_I(log(U_0mean));
end


load ../data/DataforMatlab/Compensated_base.mat
figure('name','Figure 10','NumberTitle','off')
plot(1974:2017,bt_U_med,'LineWidth',2);hold on 
plot(1974:2017,bt_U_med_base,'--','LineWidth',2); 
plot(1974:2017,bt_I_med,'LineWidth',2);
legend('Compensated $B_X$ (unobserved prices)','Compensated $B_X$ (observed prices)','Uncompensated $B_X$','Location','southeast','interpreter','latex')
%matlab2tikz('../fig/Hicksian_Mar_Median_OLS.tex');
figure('name','Figure 9','NumberTitle','off')
plot(1974:2017,bt_U_90,'LineWidth',2);hold on 
plot(1974:2017,bt_U_med,'LineWidth',2);
plot(1974:2017,bt_U_10,'LineWidth',2);
legend('90th percentile','50th percentile','10th percentile','Location','southwest','interpreter','latex')
%matlab2tikz('../fig/Hicksian_OLS.tex');


    function F_X = interpNaN(X,Y)
        indexNaN = isnan(X);
        X(indexNaN) = [];% Eliminate NaN
        Y(indexNaN) = [];% Eliminate NaN
        [X, indexUN] = unique(X); 
        F_X = griddedInterpolant(X,Y(indexUN));
    end
